###########################################################################################################
# Replication for county-level analysis presented in Table A14 (interdenominational marriages) 
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16) 
###########################################################################################################

rm(list = ls())

# Load packages
load <- c("openxlsx", "sandwich", "lmtest", "stargazer")
lapply(load, library, character.only = TRUE)


#######################################################################################
## Table A.14. Religious diversity and the prevalence of Catholic-Protestant marriages
#######################################################################################

#Read in new data at stadtkreis level
dat<-read.xlsx("county_data/refugees_intermarriage_kreise.xlsx")

#Computing key variables:
#We use only shares of Catholics and Protestants for fractionalization index because we only have data on Catholic-Protestant intermarriage rates. 
dat$Frac1939<-1-((dat$Katholische1939/(dat$Katholische1939+dat$Evangelische1939))^2+(dat$Evangelische1939/(dat$Katholische1939+dat$Evangelische1939))^2)
dat$Frac1950<-1-((dat$Kath1950/(dat$Kath1950+dat$Evang1950))^2+(dat$Evang1950/(dat$Kath1950+dat$Evang1950))^2)

dat$DeltaFrac<-dat$Frac1950-dat$Frac1939
dat$lnPop1939<-log(dat$Population1939)

dat$Popdens1939<-dat$Population1939/dat$AREA50

dat$ShareAgr1939<-dat$AgricultureInsgesamt1939/dat$Erwerbspersonen_Total1939
dat$ShareDestroyed<-dat$KriegsschaedenNWohnungen/dat$NormalWohnungen 
dat$lnDistEastBorder<-log(dat$distance_to_borderkm)

dat$ShareExp1950<-dat$Heimatvertriebene1950/dat$Pop1950


reg1<-lm(ShareMixedMarriages.1955_1958~DeltaFrac+Frac1939, data=dat)
summary(reg1) 

reg2<-lm(ShareMixedMarriages.1955_1958~DeltaFrac+Frac1939+ShareExp1950+pfemale1939+ShareAgr1939+ShareDestroyed+lnDistEastBorder+Popdens1939+lnPop1939, data=dat)
summary(reg2) 

reg3<-lm(ShareMixedMarriages.1955_1958~DeltaFrac+Frac1939+ShareExp1950+pfemale1939+ShareAgr1939+ShareDestroyed+lnDistEastBorder+Popdens1939+lnPop1939+Land, data=dat)
summary(reg3) 

reg3b<-lm(ShareMixedMarriages.1955_1958~DeltaFrac+Frac1939+pfemale1939+ShareAgr1939+ShareDestroyed+lnDistEastBorder+Popdens1939+lnPop1939+Land, data=dat)
summary(reg3b) 

reg4a<-lm(ShareMixedMarriages.1955_1958~DeltaFrac+Frac1939+ShareExp1950+pfemale1939 + ShareAgr1939+ShareDestroyed+lnDistEastBorder+Popdens1939+lnPop1939+Land, data=subset(dat, pfemale1950>median(dat$pfemale1950)))
summary(reg4a)

reg4b<-lm(ShareMixedMarriages.1955_1958~DeltaFrac+Frac1939+ShareExp1950+pfemale1939 +ShareAgr1939+ShareDestroyed+lnDistEastBorder+Popdens1939+lnPop1939+Land, data=subset(dat, pfemale1950<median(dat$pfemale1950)))
summary(reg4b)


ses1 <- list(coeftest(reg1, vcov = vcovHC(reg1, type="HC1"))[,2], coeftest(reg2, vcov = vcovHC(reg2, type="HC1"))[,2], 
             coeftest(reg3, vcov = vcovHC(reg3, type="HC1"))[,2], coeftest(reg3b, vcov = vcovHC(reg3b, type="HC1"))[,2], 
             coeftest(reg4a, vcov = vcovHC(reg4a, type="HC1"))[,2], coeftest(reg4b, vcov = vcovHC(reg4b, type="HC1"))[,2])

pvals1 <- list(coeftest(reg1, vcov = vcovHC(reg1, type="HC1"))[,4], coeftest(reg2, vcov = vcovHC(reg2, type="HC1"))[,4],  
               coeftest(reg3, vcov = vcovHC(reg3, type="HC1"))[,4], coeftest(reg3b, vcov = vcovHC(reg3b, type="HC1"))[,4], 
               coeftest(reg4a, vcov = vcovHC(reg4a, type="HC1"))[,4], coeftest(reg4b, vcov = vcovHC(reg4b, type="HC1"))[,4])


stargazer(reg1, reg2, reg3, reg3b, reg4a, reg4b, se = ses1, p=pvals1, digits=2, 
          covariate.labels=c("Delta diversity", "Religious fractionalization (1939)", "Share expellees","Share female (1939)",  "Share in agriculture (1939)", "Share destroyed",  "Ln dist to the Eastern border", "Pop density (1939)", "Ln population (1939)"), 
          omit=c("Land", "Constant"), omit.stat=c("ser","f", "rsq"), no.space=TRUE, 
          column.labels = c("Interfaith marriage, average in 1955-58", "", "","", "Above median women", "Below median women"), 
          title="Interfaith marriage. Robust SEs in parentheses.") 


